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We study the three-dimensional bosonic t-J model, i.e., the t-J model of "bosonic electrons", at 
finite temperatures. This model describes the s — ^ Heisenberg spin model with the anisotropic 
exchange coupling J± = —aJz and doped bosonic holes, which is an effective system of the Bose- 
Hubbard model with strong repulsions. The bosonic "electron" operator Bra at the site r with a 
two-component (pseudo-)spin a{— 1,2) is treated as a hard-core boson operator, and represented 
by a composite of two slave particles; a "spinon" described by a Schwinger boson (CP^ boson) 
Zra and a "holon" described by a hard-core-boson field (j)r as Bra — (j^lzra- By means of Monte 
Carlo simulations, we study its finite-temperature phase structure including the a dependence, the 
possible phenomena like appearance of checkerboard long-range order, super-counterflow, superfluid, 
and phase separation, etc. The obtained results may be taken as predictions about experiments of 
two-component cold bosonic atoms in the cubic optical lattice. 

PACS numbers: 67.85.Hj, 75.10.-b, 03.75.Nt 



I. INTRODUCTION 

Cold-atomic systems are one of the most intensively 
studied topics not only in atomic physics but also in 
condensed matter physics in these days. In particular, 
cold atoms put on an optical lattice (OL) may be used 
as a "simulator" to study certain canonical models of 
strongly-correlated electron systems^. For systems in the 
OL, interactions between atoms, dimensionality of sys- 
tem, etc. are highly controllable, and effects of impurities 
are strongly suppressed. Therefore, cold atomic systems 
in the OL are sometimes called final simulators. Among 
them, systems of double-species atoms are quite interest- 
ing from the view point of the high-temperature (T) su- 
perconductivity (SC). Investigation of these atomic sys- 
tems is expected to give an important insight into mecha- 
nism of SC in systems in which only repulsive interactions 
between particles exist. 

In this paper, we shall study the t-J model of hard- 
core bosons in the cubic lattice. There are (at least) 
two versions of the bosonic t-J model^. In the previous 
paper—, we studied one version that is a bosonic coun- 
terpart of the original fermionic t-J model and respects 
the SU(2) spin symmetry. On the other hand, in the 
present paper, we shall consider the second version that is 
an effective model of the two-band Bose-Hubbard model 
with strong repulsions and the total filling factor not ex- 
ceeding unity^. Relation and differences between these 
two versions of the bosonic t-J model were explained 
in the previous paper—. Obtained results for the sec- 
ond version of the t-J model in the present paper can 
be regarded as predictions about the system of bosonic 
atoms of two-species. Related Hubbard model at com- 
mensurate fillings has been studied in e.g., Ref.^ by the 
mean-field-theory (MFT) type approximation, and its 
one-dimensional counterpart by the Tomonaga-Luttinger 
liquid theory in Ref.^. In the present paper, we shall 



study the three-dimensional (3D) t-J model at fractional 
fillings mostly by means of the Monte- Carlo (MC) sim- 
ulations. Results are compared with the ones obtained 
previously. 

The paper is organized as follows. In Sect. 2 we ex- 
plain the model and its basic properties. We also study 
it by MFT briefly. In Sect. 3 we present the results of MC 
simulations. Section 4 is devoted for conclusions and dis- 
cussions. 



II. MODEL 

A. The t-J model 

The t-J model is derived from the Bose-Hubbard 
model^ whose Hamiltonian is given as 
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^Hub = ^ [ - + ^r+i^r + h.C.) + UnarTlbr 

r i=l 
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(2.1) 



A=a,6 



X=a,b 



where r denotes site of the cubic lattice, i{= 1,2,3) is 
the unit vector in the i-th direction (it also denotes the 
direction index), and and br are boson annihilation 
operators, nx is the number operator of the boson A, 
and therefore U and V are inter- and intra-species in- 
teractions, respectively. This i^Hub describes the sys- 
tem of two-species of cold bosonic atoms in a cubic OL. 
From Eq. ()2.ip , it is obvious that a and b atoms have the 
same hopping amplitude and the same density pa = pb in 
the present system. Recently studied ^^Rb -^^Rb atomic 
system^ is a typical example described by this Hamilto- 
nian. 

Some related models to i^Hub in Eq. (|2.ip have been 
studied so far. In the present paper, we consider the 
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specific case such that t ^U^V and the total number of 
bosons at each site is not exceeding unity (0 < Uar-^^br ^ 
1). It is obvious that the model in the above parame- 
ter region is closely related with the high-Tc materials 
and therefore it is expected that study on it gives rise to 
an important insight into the physical properties of the 
high-Tc materials. It should be remarked that at present 
the properties of the fermionic t-J model are poorly un- 
derstood in spite of the quite intensive studies on it for 
more than two decades. This fact mainly stems from the 
difficulties of numerical study on fermionic systems. 

The effective Hamiltonian in the large on-site repulsion 
limit can be derived by the standard methods of expan- 
sion in powers of t/F^ as follows; 



It is easily verified that Eq. ()2.5p is satisfied by Eq. (j2.7p . 

The hard-core boson itself can be expressed in terms 
of another CP^ boson Wrf as follows^, 



WI2W, 



• 1, ( ^r/^r/ - l)|phys)^ = 0. (2.8) 



From Eq. ()2.8p . it is obvious that lO)^ = wl2\0)w and 
^J|O)0 = ^^^Jl|0)^(;, where \0)(f){\0)uj) is the empty state of 
<j)r{wr). It is straightforward to verify that satisfies 
the mixed commutation relations of hard-core bosons^. 
Then the Hamiltonian Htj can be expressed in terms of 
the two sets of CP^ bosons Zra and Wrf. The partition 
function Z at finite T is given by the path-integral as 
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Htj = -tYl + ^r+z^r + h.c.) + J. ^ S',^,S', Z = J [DwDz] exp[- dT{zz 



(2.2) 



where the SU(2) pseudo-spin operator is given as 5*^ = 
^B^aBr with Br = {ar^brY ((J is the Pauli spin matri- 
ces). The exchange couplings are 



2r 

, J± = -Jf^--- ^ = -^Jz. (2.3) 

and Jlc is the chemical potential of holes. In the following 
discussion, we shall treat t, and Jx, hence a, as free 
parameters, and study the t-J model ()2.2p . After ob- 
taining the critical couplings etc, we shall return to the 
expression ( |27 



Physical-state condition: Double-CP^ 
representation 



^ww^ Htj{w, w, z, z))] , 
(2.9) 

where r is the imaginary time, /3 = 1/(/cbT) and 
Htj{w^w^z^z) is obtained from Htj in (|2.2p by substi- 
tuting the double-CP^ representation for and br. In 
the present numerical study, we ignore the r-dependence 
of z{r) and w{r) and consider the following system. 



Z' = J [DwDz] exp[-PHtj{iv,w,z,z)], 



(2.10) 



where z and w represent the zero- modes of z{r) and 
w{r). This approximation is justified when we consider 
system at sufficiently high temperature. However, we ex- 
pect that the system (|2.1Qp has at least qualitatively the 
same phase structure to that of (|2.9p for T > 0. As we 
discussed in RefJ>, the nonzero- modes of z{r) and w{r) 
renormalize Hfj and this renormalization tends to order 
the system. Therefore it is expected that ordered phase 
found in the system ()2.1Qp also exists in the system (|2.9p . 
This expectation was actually verified in some systems^. 



In the system of Htj in Eq. (j2.2p . a physical state at 
each site r is expanded by three orthogonal basis state 
vectors {|0), \a) = at|0), \b) = 4|0)} (a^|0) = ^0) = 0). 
In order to express this constrained Hilbert space faith- 
fully, we use the following slave-particle representation, 

ar = ^Icrl, br = (l)lcr2, (2.4) 
+ cl^Crl + cl2Cr2 - |phys) = 0, (2.5) 

where is a hard-core boson and Cra (c" = 1,2) is an 
ordinary boson. The three basis states are expressed in 
terms of Cra and \Q) {cra\^) = 0) as 

\0)^4>l\Q), al\0)^cl,\n), bl\0)^cUn). (2.6) 



In order to express the local constraint (j2.5p in more 
convenient way, we introduce a CP^ boson (Schwinger 

boson) Zra^ 

Cra = (1 - ^J0r)^ra, {cT = 1, 2) 

( ^ 4,z,,-l)|phys), = 0. (2.7) 

cr=l,2 



C. Mean- field theory 

Before going into the details of numerical study of 
Eq. (|2.1Qp . it is useful to investigate the ground-state 
properties of the model by the MET. We use a varia- 
tional wave function of bosons and br that has a site- 
factorized form, 

I*) = n I ( + f^O + |] 

r 

(2.11) 

Here we assume the sublattice symmetry and put Or = 0^ 
Xr = Xa(b) ^ A(B)-sublattice]. Then the mean- field 
energy £tj is given as 

^tJ ^ . 2/)/^ • Xa . Xb , Xa Xb\ 

— — = — sm u sm — sm h cos — cos — 

Nl 2 \ 2 2 2 2/ 



Jz 



J± 



— sm -cosxacosxb - ^sm -smxAsmxB 



-/ic COS 



(2.12) 
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where is the number of Hnks in the system. By min- 
imizing £tj^ we can obtain the MF phase diagram. The 
case in which the fihing is unity, Tij-a + rirh — 1, corre- 
sponds to ^ = TT in (|2.1ip . and 



f(x) 



Jz 

— COS XA COS xb 



^sinxAsinxB. (2.13) 



The lowest-energy state there is easily obtained as 
For J^>J^, (Xa,Xb) = (0,7r), or (xa, Xb) = (tt, 0), 
For Jx> J. >-J±, (xa,Xb) = (|,|). (2.14) 

Then for > J± , the lowest-energy state is the checker- 
board state of particle a and 6, whereas for J ^ > Jz 
the state of super-counter-flow (SCF) (albr) 7^ is real- 
ized as expected. The checkerboard state corresponds to 
an antiferromagnetic (AF) state, whereas the SCF cor- 
responds to a XY-ferromagnetic state in the magnetism 
terminology. 

Doping holes shifts to < tt. From Eq.(j232]), the 
lowest-energy state is obtained as 



For J± > Jz.XA 



XB = -, 



For Jz > J±,Xa(xb) = e, Xb(xa) = tt - e, 

-(1 - Ha- rib). {2.15} 



Jz — Jl_ 



■ cos^ - 



Jz — Jl. 



Thus, for J± > Jz Bose-Einstein condensation (BEC) 
of both a and b atoms occurs in addition to the SCF. 
On the other hand, for Jz > J±, superfluidity (SF) with 
checkerboard symmetry, so called super solid (SS), ap- 
pears for an arbitrary small but finite value of t and the 
hole densityiS. However, this result by the MFT is not 
reliable even for the present three-dimensional system be- 
cause fluctuations of the relative phases of |0), al\0) and 
bl\0) have been ignored in the MFT. 
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FIG. 2. (Left) C at t = as a function of ci for a — —0.5. 
It indicates the existence of a second-order phase transition 
(PM AF) at ci 1.51. (Right) The scaling function /(x) 
for C obtained by finite-size scaling (FSS) hypothesis^^ The 
critical exponent of the correlation length is v — 0.70. 



In the following section, we shall study the model by 
means of the MC simulations. The numerical study gives 
reliable result for the phase structure of the model and 
also details of its critical behavior. 



III. RESULTS OF MC SIMULATIONS 

A. Case of t = 

Let us turn to the numerical study We conisder the 
cubic lattice with its linear size L up to 20 and impose 
the periodic boundary condition. In order to find phase 
transition lines, we calculate the internal energy U and 
the specific heat C defined as 

U = j^{Htj), C = ^{{Hu - Uf), N = L\ (3.1) 

Furthermore we calculate various correlation functions to 
identify each observed phase. 

It is convenient to use the following parameterization, 

a = ci = pjz, C3 = pt. (3.2) 

J z 

We first consider the case of vanishing hole hopping, i.e., 
t = 0. Phase diagram was obtained for various values of 
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FIG. 1. Phase structure for t = and /ic = in the a — ci 
plane, where a — —J±/Jz and ci = I^Jz- There are three 
phases and their physical meaning is explained in the text. 
Typical value of p = (aja^) = {b^br) in each phase is also 
shown. 



FIG. 3. (Left) C as a function of ci for a = —1.5. Result 
indicates the existence of second-order phase transition (PM 
FM) at ci 1.015. (Right) Scaling function f{x). Critical 
exponent is obtained as u = 1.0. 
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FIG. 4. U and C at t = as a function of a for ci = 3. Phase 
transition takes place at a ^ —1.0. 



the chemical potential. The result for /ic = is shown 
in the a — ci plane in Figdl In the following, we shall 
mostly show results for /i^ = 0. Some of calculations of U 
and C, which were used to determine the phase bound- 
aries in FigdJ are shown in FigsISl El and |H In high-T 
region that corresponds to small ci, the system exists in 
the phase without any long-range order (LRO), which we 
call paramagnetic (PM) phase. As ci is increased, phase 
transition to ordered states takes place. For a > — 1, AF 
state with checkerboard configuration of atoms a and b 
appears as a result of strong intra-repulsion. On the other 
hand, for a < —1, the XY- ferromagnetic state appears 
at low T as a result of strong inter-repulsion. In the 
XY- ferromagnetic state, the nonvanishing condensation 
of (albr) takes place (SCF). The line a = —1, corre- 
sponding to V = 2/7, is very specific as the symmetry 
of pseudo-spin degrees of freedom is enhanced to SU(2) 
along this line, otherwise the symmetry is U(l)xZ2, i.e., 
a global {S!^ — S^) rotation and —S^ reflection. 

In the study of ferroelectric materials, the corresponding 
line is called morphotropic phase boundary (MPB), and 
it plays an important role^^. Our calculation in FigH] 
shows that the phase transition at a = — 1 looks neither 
first order nor second order. The origin of this peculiar 
behavior of U and C across the MPB is the enhancement 
of the symmetry at a = —1 as explained. Turning on the 
hole- hopping t reduces the symmetry at a = — 1 down to 
U(l)xZ2, and as a result, the phase transition becomes 
second-order. We have studied case of several other val- 
ues of /ic, and obtained a similar phase diagram to that 
in FigdJ 

The above interpretation of the phase structure is sup- 
ported by calculating the pseudo-spin correlation func- 
tions, 

^z{r) = y^(^ro^r+ro)7 

C^yi^) = jlY. E (SlS-J+r,), (3.3) 

which are used for identification of each phase (see later 
discussion and FigE]). 

To understand the properties of each phase in an intu- 
itive manner, it is helpful to examine typical configura- 
tions of variables. In FigEl we present snapshots of three 
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FIG. 5. Snapshots of three densities pa, Pb, and ph of Ea. (|3.4p 
in a XY plane in three phases of Fig.l for L — 24. From the 
above, (ci,a) = (3.0,-0.5) (AF phase), (ci,a) = (3.0,-1.5) 
(FM phase), and (ci,a) = (0.5, -0.5) (PM phase). 



densities. 

Pa {alar), ph (blbr), ph (0r0r), (3.4) 

at each phase. They are consistent with our previous 
interpretation of each phase given in the explanation of 
Fig.l. In the AF phase, atoms a and b form the checker- 
board configuration. In the FM state, the both atoms 
a and b have rather homogeneous density, and the hole 
density is very low as the energy dominates over the en- 
tropy at low T. On the other hand, the PM phase has a 
lower atomic density as the entropy dominates over the 
energy at relatively high T. 



B. Superfiuid 

In this subsection, we shall consider the case of finite 
hopping amplitude t. We verified numerically that the 
global phase structure of Fig(T] remains intact for small 
t (i.e., cs). However as cs is increased, phase tran- 
sition to SF state takes place at some critical values 
C3 = C3c(ci,a). The transition from the AF phase at 
C3 < csc to the SF phase at cs > csc is of strong first 
order as U and the hole density ph in Figl6] show. We 
employed the specific update methods for the MC sim- 
ulations in order to generate pre-choice configurations 
efficiently for the first-order phase transition-ii^. Never- 
theless, the obtained U and the hole density ph exhibit 
large hysteresis loops as cs varies. 

To verify that the SF state is realized for cs > c^o 
we calculated the boson correlation function Gair) and 
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FIG. 6. U and the hole density ph as a function of cs for 
a = —0.5 (upper panels) and 0.5 (lower panels) with ci = 3.0. 
Hysteresis loop indicates the transition from AF state to SF is 
of first-order. Similar hysteresis loops are obtained for other 
values of a and ci for the AF phase. 
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G5(r), 



Ga{r) = -^^(4o^l,ro4,ro+r^ro+r), 



Gb{r) = Xl^^^o^2,ro4,ro+r^ro+r), (3.5) 



FIG. 8. (Top) Pseudo-spin correlation functions Cxy{r) and 
Cz{r) in the AF and SF states for a = 0.5 and a = 3.0. In SF 
state, a FM LRO appears in the xy-plane of the pseudo-spin 
space. (Middle) Cxy{r) and Cz{r) in the PM and SF states 
for a = —1.5 and ci = 0.5. (Bottom) Cxy{r) and Cz{r) in the 
FM and SF states for a = —1.5 and ci = 3.0. 



and if Gai^)^ Gb{r) finite as r oo, the SF is re- 
alized. The results are shown in FiglTl It is obvious 
that Ga{r) = Gb{r) in the present case, and it has a 
nonvanishing LRO for cs = 20 indicating existence of a 
finite density of SF. In Fig 15] we also show the calcula- 
tion of the pseudo-spin correlation functions Cz{r) and 
Cxy{r). The results show that the phase transition to the 
SF state accompanies a transition from the Ising-like AF 
LRO to the XY-FM LRO^^ This result is in sharp con- 
trast with the result obtained by the MFT. The present 
numerical study indicates that the SS phase predicted in 
MFT, in which the AF LRO and SF coexist, does not 
appear in the present model. As the phase transition 
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FIG. 7. Boson correlation function G{r) — Ga(r) — Gb{r) for 
ci = 3.0, a = 0.5. Finite LRO for cs = 20 indicates that the 
phase transition in Fig|6]is a SF transition. 



to the SF phase takes place at cs = /3t ^ 0(10) and 

ci[~ 0(^,^)] = 3, the critical region is located at 
t//7, t/V ~ 3/10 <C 1 in the original Hubbard model. 
Therefore the above obtained results in the t-J model 
are also applicable for the bosonic Hubbard model. 

There are two kinds of SF, one made of atom and 
the other made of br. It is interesting to see how each SF 
behaves in the hysteresis region of the first-order phase 
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FIG. 9. Snapshots of three densities pa, pb, and ph of Ea. (|3.4p 
in a XY plane for L = 24, a = 3.0, cs = 10.0, p = 0.0, a = 0.5. 
The upper ones are in the FM+SF phase with the higher 
energy U /N — —4.39 and the lower ones are in the AF phase 
with the lower energy U /N — —6.68. 
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FIG. 10. (Upper panels) (Left) Initial configuration with 
phase separation and (Right) configuration after 60 x 10^ 
sweeps at the critical point. Dark regions represent SF with 
relatively high hole density and bright regions represent AF 
solid. (Lower panels) Densities pa, ph and ph in a horizontal 
plane in configuration after 60 x 10^ sweeps. 
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FIG. 11. Transition from PM to FM+SC. Specific heat for 
ci = 0.5 and a = — 1.5 vs C3. The peak of C has a systematic 
L dependence of a typical second-order transition. 




FIG. 12. Boson correlation function G{r) = Ga{r) — Gb{r) 
for ci — 0.5, a — —1.5 (Left) and a =3, a = —1.5 (Right). 
Finite LRO for C3 = 15 indicates that the phase transition in 
Figsim and [T3l is a SF transition. 
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FIG. 13. Specific heat of phase transition from FM to 
FM+SC, ci = 3.0 and a = —1.5 exhibits rather peculiar be- 
havior, whereas specific heat of each term ci, C3 shows typical 
system-size dependence of second-order phase transition. 



transition. In Figl9]we present snapshots of typical con- 
figurations of Pa and pb for C3 = 10.0 to check the possi- 
bility that the AF solid and SF exist separately in every 
state of update. We found that on the cs-decreasing line 
of the hysteresis loop the pure FM+SF state is realized, 
whereas the pure AF state exists on the ca-increasing 
line. This indicates that in real experiments there ex- 
ists a genuine phase transition point in the middle of 
the hysteresis loop in the MC simulation and the inter- 
nal energy has a sharp discontinuity at that transition 
point. At the discontinuity point, immiscible state of the 
AF solid and SF is realized. In order to verify this ex- 
pectation, we performed MC simulation starting with a 
half-AF and half-SF configuration and searched a "gen- 
uine critical coupling" csc at which this phase-separated 
configuration is stable during MC update. For a = 0.5 
and ci = 3.0, we found csc = 11.36, see FigHOl This re- 
sult, which shows that the phase separation takes place 
in the present 3D system, is consistent with the result of 
the previous study on the system at T = 0^. 

Let us turn to the PM SF transition. In FiglTTI 
we present C for ci = 0.5 and a = —1.5. C exhibits a 



sharp peak at C3 c^i 4.6, which indicates existence of a 
second-order phase transition. We calculated the boson 
correlation function and verified that a SF appears for 
C3 > 4.6. See FigHSti. 

We also verified that a transition from the FM to 
FM+SC takes place as C3 is increased. In the critical 
region, the total specific heat C exhibits rather pecuHar 
behavior, but the "specific heat" of each term, defined 
by Ci = {{Ei — {Ei)Y)/L'^ for each term Ei in the Hamil- 
tonian, shown in Fig (13] exhibits typical behavior of the 
second-order phase transition. In Fig|T2b. we show the 
boson correlation function for Fig(T3l The existence of 
the finite LRO means that the phase transition in FigIT3l 
is again a transition to SF. 

We observed that all three phases at t = 0, i.e., PM, 
AF and FM phases, evolve into the SF state as t is in- 
creased. Then, it is quite interesting to see if there is a 
phase boundary between these SF's for sufficiently large 
C3 though all of three phases belong to the FM+SF phase. 
This problem is closely related with recent experiment 
investigating two species SF^. This experiment observed 
that (im)miscibility of two SF's depends on the inter and 
intra-interactions between atoms. In Fig JT^ we show the 
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FIG. 14. Specific heat of eacfi term as a function of a for 
ci = 3.0 and cs = 20. There exist no anomalous behavior 
at a ^ —1. Peaks in the specific heats at a ^ —3.2 mean a 
phase transition to vanishing SF, for hole density is vanish- 
ingly small and boson correlation vanishes for a < —3.2. 



specific heat of each term and particle density as a func- 
tion of a for ci = 3.0 and cs = 20. On may expect that 
there is a phase boundary separating two FM+SF phases 
at a — 1, but the result exhibits no anomalous behav- 
iors there. On the other hand, the peaks at a o:^ —3.2 
accompanies abrupt decrease of the hole density. This 
indicates that there exists a phase transition and that 
is a transition into the vanishing SF. Ga{r) and Gb{r) 
in Fig(14] support this interpretation because they have 
no RLO at a = —4.0. The phase a < —3.2 is a pure 
FM state without holes. We also studied whether phase 
transition between two SF's takes place as the value of 
ci is varied, but we found a similar result to the above as 
varying a, i.e., there exists no phase transition between 
two SF's. 



tions. We considered the system with filling factor up 
to unity, and obtained the global phase diagram in the 
grand-canonical ensemble (GCE). At vanishing hopping 
amplitude, there are three phases in the a — ci plane, 
PM, AF and FM phases. As the hopping amplitude 
is increased, all three phases evolve into SF state with 
BEG of atoms. These obtained results are globally con- 
sistent with those for the case of integer fillings obtained 
by MFT-type approximation and numerical methods^ii^. 
However, we verified that the SS state, which is predicted 
to appear by the MFT, does not exist in the present 
model in the GGE. On the other hand, we found that 
the phase separation of the AF solid and SF is realized 
at the phase transition point. 

We also studied if there exists phase boundary between 
the SF's. However there are no phase boundaries between 
them. 

Results obtained in the present paper show that the 
bosonic t-J model has a very rich phase structure. We 
studied the system in the GGE. It is quite interesting to 
study the bosonic t-J model in the canonical ensemble 
with fixed average atomic number. In particular, an in- 
homogeneous state may appear near the first-order phase 
transition point from the AF solid to the SF. This prob- 
lem is under study and results will be reported in a future 
publication. 
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